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In the context of the recently developed "equation-free" approach to computer-assisted analysis 
of complex systems, we extract the self-similar solution describing core collapse of a stellar system 
from numerical experiments. The technique allows us to side-step the core "bounce" that occurs in 
direct A^-body simulations due to the small- iV correlations that develop in the late stages of collapse, 
and hence to follow the evolution well into the self-similar regime. 
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I. INTRODUCTION 

In many areas of current research, physical models are 
available at a fine, microscopic scale, while the questions 
of most interest concern the system behavior on a coarse- 
grained, macroscopic level. An example is the gravi- 
tational A^-body problem p|. Coarse-grained approxi- 
mations exist, e. g. the orbit-averaged Fokker-Planck 
equation , which approximates the full 6N equations of 
motion by a differential operator that acts on the single- 
particle distribution function /. But the derivation of the 
Fokker-Planck equation from the full A^-body equations 
of motion requires a number of approximations. The 
"equation-free" computational framework [J- [3- Q has 
recently been proposed for the computer-assisted study 
of such problems, circumventing the derivation of explicit 
coarse-grained equations. The coarse-grained behavior is 
analyzed directly, through appropriately designed short 
computational experiments by the fine-scale models. In 
the case of problems where the macroscopic behavior is 
characterized by scale invariance, dynamic renormaliza- 
tion 0, S 13 1 combined with equation- free computation 
and a template-based approach 0, can in effect convert 
the self-similar problem into a stationary one, by work- 
ing in a frame of reference that expands or shrinks along 
with the macroscopic system observables. 

Here we apply equation-free dynamic renormalization 
techniques to the problem of gravitational core collapse 
[Tol Hlj . A star cluster or galaxy evolves due to grav- 
itational encounters between the stars, which drive the 
local velocity distribution toward a Maxwellian. Stars in 
the high- velocity tail of the Maxwellian escape from the 
system; the probability of escape in one relaxation time 
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is roughly 1%. Here is the mean square (3D) ve- 
locity of the stars, p is the mass density, m is the indi- 
vidual stellar mass, G is the gravitational constant, and 



In A ps ln(0.4iV) is the Coulomb logarithm @, with N 
the number of stars in the cluster. Escape occurs pri- 
marily from the high-density central region, or "core" ; if 
the density contrast between core and envelope is suffi- 
ciently great, the core exhibits the negative specific heat 
characteristic of gravitationally bound systems and 
contracts. 

Treatments of core collapse based on fluid [l3L IhH ] or 
Fokker-Planck pH , Il6l IT^| approximations to the full N- 

body equations of motion suggest that the late stages are 
self-similar, 
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with p c and r c the central density and core radius re- 
spectively, and t cc the time at which the central density 
diverges; t cc — t is roughly 330 times the relaxation time 
of Eq. I if the latter is evaluated at the center of the 
cluster. 

These predictions are in reasonable a gree ment with the 
results of direct A^-body integrations jisl llH |2(j . But 
when N is finite, the number of particles in the core 
drops to zero as collapse proceeds, causing two-body and 
higher order correlations to dominate the evolution; typ- 
ically, binary stars form which halt the collapse ("core 
bounce"). In all A^-body simulations published to date, 
this occurs before or only shortly after the core has en- 
tered the self-similar regime. The rate of binary forma- 
tion per unit volume due to three-body interactions is 
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By the time that the number of stars in the core has 
dropped to N c (N c = p c r c)i the total number of bina- 
ries formed is ~ lO 3 ^" 2 . Since a single hard binary can 
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halt the collapse, bounce occurs when N c has dropped 
to ~ 30. This limits the maximum central density that 
can be reached in an JV-body simulation to a multiple 
~ 10 35 (./V/10 4 ) 3 of the initial density (assuming the ini- 
tial state described below); for N ss 10 4 , the achievable 
density contrast is ~ 10 4 , which is also roughly where 
self-similar behavior first appears [Lsl ] . 

In this paper, we exploit coarse dynamic renormal- 
ization to compute representative self-similar solutions 
at scales "away from" the latest stages of core collapse. 
This allows us to avoid the finite- N correlations that de- 
velop at those stages and recover features of the self- 
similar behavior predicted at the large- TV limit. We de- 
fine the "macroscopic" quantity of interest to be the 
single-particle distribution function f(E), E — v 2 /2 + 
^f(r), with <3>(r) the self-consistent gravitational poten- 
tial (spherical symmetry is assumed throughout). By 
dynamically renormalizing the A-body model after each 
short "burst" of integration, we are able to maintain a 
large (~ 10 3 ) number of particles in the core even with 
modest (~ 1QK) total numbers, effectively halting binary 
formation and allowing us to accurately extract the form 
of the self-similar solution. 



II. DESCRIPTION OF THE CALCULATIONS 



via Poisson's equation. The isotropic, single-particle dis- 
tribution function f(E) corresponding to this density- 
potential pair is given uniquely by Eddington's formula 
|24j. For the purposes of generating a new Monte-Carlo 
set of positions and velocities, the relevant function is 
not f(E) but N(< v,r), the cumulative number of stars 
with velocities less than v at radius r. Using Eddington's 
formula, this can be shown to be 
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2. Integrate. The TV-body realization was advanced in 
time until the central density had increased by a factor of 
~ 5 compared with the initial model. In the nearly-self- 
similar regime, this required approximately 250 central 
relaxation times. 

3. Restrict. An estimate of p(r) was computed from the 
particle coordinates x», i = 1...N at the final time step. 
The position of the cluster center was first determined as 
in [25J and the particle radii n were defined with respect 
to this center. Each particle was then replaced by the 
kernel function 



We adopted the standard Plummer model |2S| initial 
conditions for this problem, with density, gravitational 
potential, and single-particle distribution function given 
by 
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Here and below, the gravitational constant G has been 
set to one. The iV-body algorithm was an adaptation 
of S. J. Aarseth's NB0DY1 code to the GRAPE-6 
special-purpose computer and included a fourth-order 
particle advancement scheme with individual, adaptive, 
block time steps. We first used this code to carry out 
a direct integration until core bounce of a Plummer 
model with N = 16384 particles; the results (e.g. evolu- 
tion of the central density, time of core bounce) agreed 
well with published studies using different TV-body codes 
[IB Ell H3 • We used the same number of particles in the 
calculations described below. 

Our coarse renormalization procedure consisted of 
short bursts of simulation in a "lift-simulate-restrict- 
rescale-truncate" procedure, repeated in a loop until the 
asymptotic form of the self-similar solution emerged. In 
detail, the steps were: 

1. Lift: Given a smooth representation p(r) of the 
density profile, a set of Monte-Carlo positions and ve- 
locities was generated, as follows. First an estimate of 
the gravitational potential \&(r) was computed from p(r) 
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which is an average over the sphere of radius r; of the 
3D Gaussian kernel of width The density estimate, 
P( r ) — J2iLih^ 3 K(r,ri,hi), was then computed on a 
logarithmic radial grid. The kernel width hi associated 
with the ith particle was determined by first construct- 
ing a pilot estimate of the density via a nearest-neighbor 
scheme, then allowing the Ju to vary as the inverse square 
root of this pilot density [26j . 

4. Rescale: The density estimate was rescaled, p(r) — > 
Ap(Br), according to an algorithm that left the core 
properties unchanged in the self-similar regime. The ver- 
tical normalization A was fixed by comparing the values 
of the density at r = at the start and end of the inte- 
gration interval. The radial scale factor B was adjusted 
such that the two density estimates had the same value 
at the radius where the density was 1/50 of its central 
value. 

5. Truncate: In the late stages of core collapse, the 
system develops ap~ r~ Q , a ~ 2 envelope around the 
shrinking core. In direct iV-body integrations, this en- 
velope grows to contain most of the total mass, and the 
number of stars in the core drops to zero, leading to the 
small- A effects discussed above. We avoided this by trun- 
cating the rescaled density beyond a radius r± w 50r c 
using a function that falls smoothly to zero at a radius 
r 2 w 10ri. 

6. Step 1 was then repeated in a loop. 
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FIG. 1: Evolution of the rescaled density. Dashed line is the 
initial state and points are the self-similar p*(r) derived in 
Ref. Inset shows a = — dlogp/d logr at the end of the 

final integration; dashed line is a — —2.23. 
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FIG. 2: Evolution of the dimensionless collapse rate param- 
eter £ = (p c / p c )/t r . Each time interval corresponds to one 
"burst" of integration, after which the model was rescaled as 
described in the text. Dashed line shows the asymptotic (self- 
similar) value £ = 3.6 x 10 -3 as computed via the isotropic 
orbit-averaged Fokker-Planck equation ^| 0] . 



Generation of the new particle coordinates from the 
single-particle f(E) in step 1 has the effect of removing 
any binaries that may have formed in the previous inte- 
gration interval, and of resetting to zero any anisotropics 
that developed in the velocities. Since the number of 
particles in the core after each rescaling was ~ 10 3 , the 
chance that even a single binary would form before the 
next rescaling was negligible. Velocity anisotropy was 
evaluated at the end of each interval and found to be very 
small, and restricted to the largest radii. The effects of 
binary formation and anisotropy growth could be reduced 
still more by increasing N and by reducing the lengths of 
the integration intervals. Newton-type fixed point algo- 
rithms that converge on the self-similar solution are also 
possible. 



III. RESULTS 

Fig. 1 shows the density profile at the end of each 
rescaling step, compared with the self-similar p*(r) com- 
puted from the orbit-averaged isotropic Fokker-Planck 
equation [l^, which has /?+ oc r -2,23 at large radii. The 
renormalized density is quite close to the self-similar so- 
lution after 4-5 iterations. 

In the self-similar regime, the finite changes in cen- 
tral density and core radius during one integration in- 
terval satisfy log(p c j / p Cji ) / logjrv //r C)i ) = j/6; in the 
Fokker-Planck approximation |l7|. this ratio is equal 
to the asymptotic density slope, j/6 = a « —2.23. 
We computed this ratio from the rescaling factors 
[A, B) in each integration interval and found it to 
be (-2.09, -2.15, -2.23, -2.24, -2.23, -2.23), consistent 



with the Fokker-Planck prediction at late times, and also 
consistent with our computed value of a (Fig. 1). 
The dimensionless collapse rate is 

^t r (0)-^ (7) 
Pc at 

with t r (0) the value of the relaxation time (Eq. 1) eval- 
uated at r = 0; in the self-similar regime, £ should reach 
a constant value. Fig. 2 shows the evolution of £, com- 
puted by fitting a smoothing spline to p c (t) during each of 
the integration intervals. Isotropic Fokker-Planck treat- 
ments llo . E3 find £ w 3.6 x 10 3 in the self-similar 
regime, consistent within the noise with Fig. 2. 

Another way to measure the degree of evolution 
achieved via the rescaled integrations is shown in Fig. 
3. In this plot, the cumulative effect of the rescalings on 
p c and r c is included; the central density has increased 
by a cumulative factor of ~ 10 4 5 by the end of the fifth 
rescaling, and the core radius has decreased by a factor 
~ 10 2 . The p c (r c ) relation from the iV-body integra- 
tion without rescaling is also shown. In the un-rcscaled 
integration, the central density peaks due to binary for- 
mation at a value ~ 10 times lower than in the rescaled 
runs. 

In principle, one can extract the similarity exponents 
7 and 6 of Eq. 2 (not just their ratio) from the numerical 
models. For instance, if t\ and ti are two distinct times 
in the self-similar regime, then 
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FIG. 3: Central density vs core radius in the rescaled 
runs (points) and in a direct JV-body integration (blue line). 
Dashed line shows p c oc r~ 2 ' 23 . In the direct JV-body sim- 
ulation without rescaling, the central density reaches a peak 
value then decreases ( "core bounce" ) due to formation of bi- 
nary stars. 

and similarly for S. Another approach Q is via numerical 
calculation of the scaling exponents that characterize the 
effective differential operator. A third approach is sim- 
ply to fit Eqs. 1 to the time-dependent density and core 
radius in the self-similar regime. We had limited suc- 
cess with each of these methods due to noise associated 
with the modest particle numbers. The noise could be 
reduced by averaging the results from different random 
realizations of the same initial conditions, or by increas- 
ing N. 



namical system to a stationary state, in a scaled reference 
frame where the self similarity has been "factored out." 
We have demonstrated the usefulness of the method for 
studying gravitational core collapse. In the limit of large 
particle numbers, core collapse is characterized by a cen- 
tral density that increases without limit, due to the com- 
bined effects of heat transfer and the negative specific 
heat of the core. In numerical simulations with finite 
N, collapse is halted when the number of particles in 
the shrinking core drops to a small enough value that 
binaries can form. By carrying out short bursts of in- 
tegration of appropriately rescaled initial conditions, we 
showed that it is possible to avoid these small- N effects 
and to follow collapse well into the self-similar regime. 
Furthermore we achieved this with quite modest particle 
numbers, N 10 4 . In our approach, the degree to which 
core collapse can be followed is essentially independent 
of the number of particles used, while in direct A-body 
simulations, the time to core bounce is determined by N. 

We chose the macroscopic function of interest to be 
the isotropic, single-particle distribution function f(E, t). 
The same approximation is commonly made in Fokker- 
Planck and fluid treatements of core collapse |l3lll6lll7| . 
A larger particle number would allow us to relax the 
assumption of isotropy and extract f(E,L,t) with L 
the orbital angular momentum per unit mass. Results 
so obtained could be compared with solutions to the 
anisotropic Fokker-Planck HH and fluid Q equa- 
tions. 
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